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A generic model of stochastic autocatalytic dynamics with 
many degrees of freedom Wi, i = 1, . . . ,N is studied using 
computer simulations. The time evolution of the Wi's com- 
bines a random multiplicative dynamics Wi(t + 1) = Xwi(t) 
at the individual level with a global coupling through a con- 
straint which does not allow the Wi's to fall below a lower 
cutoff given by c-w, where w is their momentary average and 
< c < 1 is a constant. The dynamic variables Wi are found 
to exhibit a power-law distribution of the form p(w) ~ ui _1_Q1 . 
The exponent a(c, N) is quite insensitive to the distribution 
n(A) of the random factor A, but it is non-universal, and 
increases monotonically as a function of c. The "thermody- 
namic" limit N — > oo and the limit of decoupled free multi- 
plicative random walks c — ► do not commute: a(0,N) — 
for any finite N while a(c, oo) > 1 (which is the common range 
in empirical systems) for any positive c. The time evolution 
of w(t) exhibits intermittent fluctuations parametrized by a 
(truncated) Levy-stable distribution L a (r) with the same in- 
dex a. This non-trivial relation between the distribution of 
the Wi 's at a given time and the temporal fluctuations of their 
average is examined and its relevance to empirical systems is 
discussed. 



I. INTRODUCTION 

The origins of power-law distributions as well as their 
conceptual implications have been an active topic of re- 
search in recent years. Power laws are intrinsically re- 
lated to the emergence of macroscopic features which are 
scale invariant within some bounds and distinct from the 
microscopic elementary degrees of freedom. Often, these 
features are insensitive to the details of the microscopic 
structures. Well known examples of power law distri- 
butions include the energy distribution between scales in 
turbulence |l[] , the distribution of earthquake magnitudes 
Q , the diameter distribution of craters and asteroids || , 
the distribution of city populations the distribu- 

tions of income and of wealth [@-{L2), the size-distribution 
of business firms 13 and the distribution of the fre- 
quency of appearance of words in texts . The fact that 
multiplicative dynamics tends to generate power-law dis- 
tributions was intuitively invoked long ago |I^ , [l5| -|l7| but 
the limitations in computer simulation power kept the 
models under the constraints imposed by the applicabil- 
ity of analytical treatment. More recently, a broader class 
of models has been studied combining computer simu- 



lations with theoretical analysis within the Microscopic 
Representation paradigm proposed in Ref. [i8fl . In par- 
ticular, it was shown p9|-pi[ that power laws appear in a 
variety of dynamical processes and are maintained even 
under highly non-stationary conditions. 

In this paper we consider a generic model of stochas- 
tic dynamics with many degrees of freedom Wi(t), i — 
1,...,N. The time evolution of the w^s is described 
by an asynchronous update mechanism in which at each 
time step one variable is chosen randomly and is multi- 
plied by a factor A taken from a predefined distribution. 
In addition, there is a global coupling constraint which 
does not allow the Wi 's to fall below the lower cutoff given 
by c • w, where w is the momentary average of the Wi's 
and < c < 1 is a constant. The dynamic variables 
Wi are found to exhibit a power-law distribution of the 
form p(w) ~ w~ l ~ a . The exponent a is found to be in- 
sensitive to the distribution 11(A) of the random factor 
A. However, a is non-universal, and increases monoton- 
ically as a function of c. In the limit c = (where the 
Wi's become decoupled) a — for any finite N. How- 
ever, in the "thermodynamic" limit N — oo, a > 1 for 
any positive c. Thus the two limits do not commute. 
This is important for applications since typically in em- 
pirical systems a > 1 |l-|l4]| unlike the case of the free 
multiplicative random walk which predicts a log-normal 
distribution J2^,^3| corresponding to a — Q . 

The time evolution of w(t) exhibits intermittent fluc- 
tuations parametrized by a truncated Levy-stable distri- 
bution with the same index a. This intricate relation 
between the distribution of the w^'s at a given time and 
the temporal fluctuations of their average is examined 
and its relevance to empirical systems is discussed. Our 
model indicates that in certain cases the scaling exponent 
may be insensitive to the distribution of the multiplica- 
tive (random) factor A and depends only on the "lower 
bound" features which control the smallest values of the 
elementary variables. The relation between the limiting 
conditions and the power law exponent is to be applied 
in each particular case and it constitutes a strong instru- 
ment in identifying and validating the relevant degrees of 
freedom responsible for the emergence of scaling. 

The present paper proposes to consolidate by numer- 
ical simulations the control one has on a specific model 
and help in this way its further application to additional 
systems. The paper is organized as follows. In Sec. II we 
present the model. Simulations and results are reported 
in Sec. Ill, followed by a discussion in Sec. IV and a 
summary in Sec. V. 
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II. THE MODEL 



A. Formal Definition 



The model [19 20] describes the evolution in discrete 
time of N dynamic variables Wi(t), i = 1, . . . , N. At each 
time step t, an integer i is chosen randomly in the range 
1 < i < N, which is the index of the dynamic variable Wi 
to be updated at that time step. A random multiplicative 
factor X(t) is then drawn from a given distribution n(A), 
which is independent of i and t and satisfies J. II(A)ciA = 
1. This can be, for example, a uniform distribution in 
the range X min < X < X rnax , where X min and X max are 
predefined limits. The system is then updated according 
to the following stochastic time evolution equation 



Wi(t+1) = \(t)Wi(t) 

Wj(t+l) = Wj {t), j 1 V: j • ;. 



(1) 



This is an asynchronous update mechanism. The average 
value of the system components at time t is given by 



1 N 



(2) 



The term on the right hand side of Eq. (|l|) describes 
the effect of auto-catalysis at the individual level. In 
addition to the update rule of Eq. the value of the 
updated variable Wi(t + 1) is constrained to be larger or 
equal to some lower bound which is proportional to the 
momentary average value of the Wi 's according to 



Wi(t + 1) > c • w(t) 



(3) 



where < c < 1 is a constant factor. This constraint is 
imposed immediately after step (pi) by setting 



Wi(t + 1) — > max{wi(t + 1), c ■ w(t)}, 



(4) 



where w(t), evaluated just before the application of 
Eq. (Q), is used. This constraint describes the effect of 
auto-catalysis at the community level. 



B. Main Features 

Our model is characterized by a fixed (conserved) num- 
ber of dynamic variables N, while the sum of their val- 
ues is not conserved. The conservation of the number of 
dynamic variables, which is enforced through the lower 
cutoff constraint is essential since otherwise the system 
dwindles over time. The non-conservation of the sum 
of the values of the dynamic variables is important as 
well. It allows to perform the multiplicative updating on 
a single variable at a time with no explicit binary interac- 
tions since a gain in uij does not require a corresponding 
immediate loss by other Wj's. In fact, the interactions 
between the dynamic variables are implied only in the 



step of Eq. (||) in which the lower cutoff is imposed. The 
dynamic rule (jl|) can be described by a master equation 
for the probability distribution p(w) of the form 



p(w,t + 1) -p(w,t) = — 



Tl(X)p(w/X,t)dX-p(w,t) 



(5) 



where the 1/N factor takes into account the fact that 
only one of the w^s is updated in each time step. This 
description applies for the bulk of the distribution of the 
uii's but not in the vicinity of the lower cutoff where the 
step of Eq. (|[) which is not taken into account by Eq. (|^) 
may be dominant. 

For the following analysis it is convenient to normalize 
the Wj 's according to 



Wj(t) -> Wj(t)/w(t), j = l, 



, N. 



As a result, the new average w(t) is normalized to 



w(t) 



N 



wp(w, t)dw = 1, 



(6) 



(7) 



while ^iWiit) — Nw — N. Performing this normaliza- 
tion step after each iteration removes the non station- 
ary part of the distribution and amounts statistically 
to an overall multiplicative factor. This (time depen- 
dent) factor which represents a global inflation rate can 
be recorded at each step. It is convenient to represent 
the dynamics (^|) on the logarithmic scale. In terms of 
the new variables 



W t = Inwi, 



(8) 



Eq. (|I|) defines a random walk with steps of random size 
In A: 



Wi(t+ 1) =Wi(t) + In A. 



(9) 



The corresponding probability distribution P(W) be- 
comes 

P(W) =e w p(e w ). (10) 
In terms of P and W, the master equation (||) becomes: 



P(W,t + l) - P(W,t) = 
1 

N 



J U(X)P(W -In X,t)dX- P(W,t) 



(11) 



The asymptotic stationary solution, is found to be fl9[l 



P(W) 



-aW 



(12) 



In terms of the original variable Wi , we get according to 
Eq. ( |l(3| ) a power law distribution: 



p(w) — K ■ w 



(13) 
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The value of the exponent a is determined by the nor- 
malization condition [Eq. (0)] divided by the probability 

normalization condition J p(w,t)dw — 1 (in order to 
eliminate the constant factor K), which yields: 



TV 



(#r-(#) 



(14) 



The exponent a is given implicitly as a function of c 
and iV by Eq. (|l4|). We identify two regimes within < 
c < 1 in which Eq. ( [Ti] ) can be simplified and a can be 
obtained explicitly. For a given TV and values of c in the 
range 1/ In TV <C c < 1 one obtains a > 1 as well as 
(c/TV) Q -C c/TV -C 1. Consequently, in this range, one 
can neglect the (c/N) a terms in Eq. (S3) to obtain to a 
good approximation 



TV = 



a — 1 



-1 



-(#) 



which gives the explicit, TV-independent solution 

1 



1-C 



(15) 



(16) 



This relation is exact in the "thermodynamic" limit TV = 
oo. The relation ([l6]) has two remarkable properties: (a) 
it does not depend on the distribution 11(A); (b) it gives 
rise to a values in the experimentally realistic range a > 
1. 

For finite TV and values of c lower than 1/ In TV the 
approximation Eq. ( |l6| ) breaks down and values a < 1 
become possible. However, for any finite N, another ap- 
proximation holds in the range c <C X/N < X. In this 
range (c/N) <C (c/N) a <C 1 and therefore one can ne- 
glect (c/TV) Q in the numerator of Eq. ( |l4| ) and c/TV in 
the denominator to obtain: 



TV 



a — 1 



-1 



(17) 



By taking the logarithm on both sides and neglecting 
terms of order 1 we obtain 



In TV 

MnJc) ' 



(18) 



Note that even for systems in which the lower bound 
(which is due to some microscopic discretization) given 
by c, is orders of magnitude smaller than X/N, the result- 
ing a may differ significantly from the free multiplicative 
random walk result a — 0. Since c enters in the formula 
( |l8| ) for a through its logarithm, the system gives away 
information on its microscopic scale cut-off c through the 
exponent a of its macroscopic power law behavior. 

One should emphasize that in the region where a < 1 
the average w of the distribution p(w) in Eq. (^) is not 
well defined and in fact one expects in the actual runs 



very wide macroscopic fluctuations of this mean. These 
fluctuations are however never infinite because accord- 
ing to the formulae above, as one increases the size of 
the system TV, the region along the c axis where a < 1 
shrinks to 0. For 1 < a < 2 it is only the standard devia- 
tion of the distribution p(w) which is formally divergent. 
This gives rise in the actual computer simulations to wide 
fluctuations of the individual values of Wi. However, this 
divergence is kept in check too by the fact that no Wi can 
possibly exceed N-w, namely p(N-w) = 0. This amounts 
to a truncation from above of the power law Eq. (Oh . 



III. NUMERICAL SIMULATIONS AND RESULTS 

Numerical simulations of the stochastic multiplicative 
process described by Eqs. ([j]) and (Q), confirm the va- 
lidity of Eq. ( p"3| ) for a wide range of lower bounds c. It 
appears that the exponent a is largely independent of 
the shape of the probability distribution n(A). Fig. |l| 
shows the distribution of Wi, i = X, . . . , N, obtained for 
N = 1000, c = 0.3, and A uniformly distributed in the 
range 0.9 < A < 1.1. A power law distribution is found 
for a range of three decades between w m i n = 0.0003 and 
Wmax — 0.3. The slope of the best linear fit within this 
range is given by a = 1.4, in agreement with Eqs. jl^ ) 
and ([l6]). On the horizontal axis of this graph the sum of 
all Wi's is normalized to 1 and therefore w = 0.001. The 
exponent a as a function of the lower cutoff c is shown 
in Fig. |[ Numerical results are presented for N = 100 
(empty dots), 1000 (full dots) and 5000 (squares). The 
prediction of Eq. (jlj) is shown for N — 1000 (solid 
line), which is in good agreement with the numerical 
results for all values of c. The approximate expression 
Eq. ( |l6| ) is also shown (dashed line). It is observed that 
for TV = 1000 this approximation gradually starts to hold 
as c is increased beyond l/ln(1000), in agreement with 
the theoretical analysis. In general, for a given TV, a is 
monotonically increasing as a function of c, starting from 
a = (which corresponds to X/w distribution) at c = 0, 
where the w;'s are uncoupled. It is also observed that 
as TV is raised, the value of a which corresponds to a 
given c increases monotonically. As a result, the range of 
validity l/lnTV^;c<lof the approximation Eq. ( |l6| ) 
is extended and the knee adjacent to c = sharpens 
and becomes a discontinuity for TV — > oo. The range 
<C c < 1 /TV in which the approximation of Eq. ( |l8| ) is 
valid, shrinks correspondingly. 

Let us turn now to the dynamics of the system as a 
whole. The dynamics of the system involves, according 
to Eq. ([j]) , a generalized random walk with step sizes dis- 
tributed according to Eq. (|l3|). Therefore, the stochastic 
fluctuations of w(t) after r time steps: 



r(r) 



i{t + T)-w(t) 
w(t) 



(19) 



are governed p5[ by a truncated Levy distribution L a (r). 
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In Fig. [3| we show the distribution of the stochas- 
tic fluctuations r(r) for r = 50, which is given by a 
(truncated) Levy distribution L a (r). According to Ref. 
|P6| , the peak of the (truncated) Levy-stable distribution 
scales with r as 

L a (r = Q) ~ r- x ' a (20) 

where a is the index of the Levy distribution. In Fig. ^ 
we show the height of the peak P(r = 0) of Fig. || as 
a function of r. It is found that the slope of the fit in 
Fig. H is —0.71, which following the scaling relation J2C|) 
means that the index of the Levy distribution in Fig. H 
is a — — 1/(— 0.71) = 1.4. These results were obtained 
for the same parameters which gave rise to the power law 
distribution with a = 1.4 in Fig. fy. Thus, the prediction 
that the fluctuations of w in Fig. Q follow a (truncated) 
Levy-stable distribution with an index a which equals 
the exponent a of the power-law distribution in Fig. ||, 
is confirmed. 



IV. DISCUSSION 

The model considered in this paper may be relevant to 
a variety of empirical systems in the physical, biological 
and social sciences which can be described by a set of 
interacting dynamic variables which follow a stochastic 
multiplicative dynamics. Such dynamical processes may 
play a role in the formation of the mass distribution in 
the universe where clusters of galaxies accumulate and 
eventually form super-clusters. In a different context, 
the growth of cities is basically a multiplicative process 
governed by the reproduction rate of the local population 
in addition to mobility between cities. 

Enhanced diffusion processes, which can be described 
by the Levy-stable distribution have been observed in 
a variety of nonlinear dynamical systems |27],|2£| . Unlike 
the stochastic model studied here, these systems are gov- 
erned by deterministic rules. They exhibit intermittent 
chaotic motion which gives rise to enhanced diffusion. 

In population dynamics, the number of individuals in 
each specie varies stochastically from one season to the 
next with a multiplicative factor which depends on the 
local conditions. The lower bound may represent the 
minimal number of individuals required for the species 
to survive in the given environment. In this case the 
number of species may not be strictly a constant, but 
species that are wiped out may be replaced by others 
which invade their area. In this context it was found 
that the number of species of a given size often follows a 
decreasing power-law distribution as a function of their 
size (see e.g. Ref. p9fl). 

In the economic context of a stock-market system the 
dynamic variables Wi, i = 1 , . . . , N may represent the 
wealth of individual investors. In this case the dynamics 
represents the increase (or decrease) by a random factor 
X(t) of the wealth wi of the investor i between times t and 



t + 1. The lower bound may represent a minimal wealth 
required in order to participate in stock market trading. 
In a more general economic model, this lower bound may 
be related to a basket of basic publicly funded services 
which every individual receives. In another possible in- 
terpretation, the WiS represent the capitalization (total 
market value) of the firm i, which may increase (or de- 
crease) by a factor \(t) at each time step. In this case the 
lower bound may represent the minimal requirements for 
a company stock to be publicly traded. 

Studies of the distribution of wealth in the general 
population revealed a power-law behavior (see e.g. Ref. 
p0{). More recently it was shown Q that the distribu- 
tion of individual wealth of the 400 richest people in the 
United States (Forbes 400) corresponds to a power law 
with a = 1.36 [more precisely W(n) — C ■ n^ 1 /" where 
W(n) is the wealth of the n-th richest person on the list]. 
Recent analysis of stock market returns, measured over 
many years found a truncated Levy distribution L a (r) 
with the index a = 1.4 for an extended (but finite) range 
of returns r [[u) . These results indicate that the property 
observed in our model, namely that the same value of the 
index a appears both in the power law distribution and 
in the Levy-stable distribution of the fluctuations may 
be of relevance in the economic context. To further ex- 
plore this possibility it would be interesting to examine 
whether the distribution of total market values of com- 
panies in the stock market exhibits a power law behavior 
of the form ( |l3| ) with a = 1.4. 

V. SUMMARY 

We have studied a generic model of stochastic auto- 
catalytic dynamics of many degrees of freedom using 
computer simulations. The model consists of dynamic 
variables Wi, i = 1 , . . . , N which are updated randomly 
one at a time through an autocatalytic process at the 
individual level. In addition, the variables are coupled 
through a lower bound constraint which enhances the 
variables which fall below a fraction of the global average. 
The model may describe a large variety of systems such 
as stock markets and city populations. The distribution 
p(w,t) of the system components Wi turns out to fulfill 
a power law distribution of the form p(w,t) ~ w^ 1 ^". 
In the limit N — oo, c — > one obtains the case often 
encountered in nature: awl. The average w(t) exhibits 
intermittent fluctuations following a Levy-stable distri- 
bution with the same index a. This relation between 
the distribution of system components and the temporal 
fluctuations of their average may be relevant to a variety 
of empirical systems. For example, it may provide a con- 
nection between the distribution of wealth/capitalization 
in a stock market and the distribution of the index fluc- 
tuations. 
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FIG. 1. The distribution of the variables w iy i = 1, . . . , N for N = 1000 obtained from a numerical simulation of the model 
given by Eqs. (|l| and (Q) with the lower cutoff c = 0.3 and 11(A) uniformly distributed in the range 0.9 < A < 1.1. The 
distribution (presented here on a log — log scale) exhibits a power law behavior described by p(w) ~ w~ 1 ~ a , where a — 1.4. 

FIG. 2. The exponent a of the power-law distribution of the variables Wi, i — 1, . . . , N as a function of the lower cutoff c. 
The data were obtained from the simulations of the multiplicative stochastic process of Eqs. (|l|) and (Q) with N = 100 (empty 
dots), 1000 (full dots) and 5000 (squares). The theoretical prediction of Eq. ( |l4| ) is shown for N = 1000 (solid line) and is in 
excellent agreement with the numerical values for all values of c. The approximate expression of Eq. ([l6]) is also shown (dashed 
line) . 

FIG. 3. The distribution of the variations of w after r steps r(r) = [w(t + r) — w(t)]/w(t) where r = 50, for the same 
parameters as in Fig. [I]. This distribution has a Levy-stable shape with index a = 1.4. 

FIG. 4. The scaling with r of the probability that r(r) = [w(t + t) — w(t)\/w(t) is 0. The parameters of the process are the 
same as in Figs. [I] and ^. The slope of the straight line on the logarithmic scale is 0.71 which corresponds to a Levy-stable 
process with a = 1/0.71 = 1.4. 
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